function CountryCell=CaballeroEngelDecomposition_func(CountryCell,ii)

global switch_figure;

deviationVec = CountryCell.deviationVec;
hazardVec = CountryCell.("freq_any_h" + num2str(ii) + "_meanVec");
densityVec = CountryCell.densityVec;
lengthA = CountryCell.lengthA;

lengthdeviation = size(deviationVec,1);

%derivative of the hazard using left and right values
LambdaPrime(2:lengthdeviation-1,1) = (hazardVec(3:lengthdeviation,1)-hazardVec(1:lengthdeviation-2,1))./(2*(deviationVec(2,1)-deviationVec(1,1)));
LambdaPrime(1,1) = (hazardVec(2,1)-hazardVec(1,1))./((deviationVec(2,1)-deviationVec(1,1)));
LambdaPrime(lengthdeviation,1) = (hazardVec(lengthdeviation,1)-hazardVec(lengthdeviation-1,1))./((deviationVec(2,1)-deviationVec(1,1)));

%Calculate values used in Karadi-Schoenle-Wursten, 2021 section 5.1
deviation0Index = find(deviationVec==0);

negdeviationRange = 1:deviation0Index-1;
posdeviationRange = deviation0Index:lengthA-2;

F0 = sum(densityVec(negdeviationRange));
gg = densityVec(negdeviationRange)/F0;
hh = densityVec(posdeviationRange)/(1-F0);

inflMinus = sum(-deviationVec(posdeviationRange).*hazardVec(posdeviationRange).*hh);
xbarMinus = sum(deviationVec(posdeviationRange).*hh);
LambdabarMinus = sum(hazardVec(posdeviationRange).*hh);

intensiveMinus = LambdabarMinus;
LambdaPrimebarMinus = sum(LambdaPrime(posdeviationRange).*hh);
grossExtensiveMinus = xbarMinus*LambdaPrimebarMinus;
selectionMinus = sum(deviationVec(posdeviationRange).*(LambdaPrime(posdeviationRange)-LambdaPrimebarMinus).*hh);


inflPlus = sum(-deviationVec(negdeviationRange).*hazardVec(negdeviationRange).*gg);
xbarPlus = sum(deviationVec(negdeviationRange).*gg);
LambdabarPlus = sum(hazardVec(negdeviationRange).*gg);

intensivePlus = LambdabarPlus;
LambdaPrimebarPlus = sum(LambdaPrime(negdeviationRange).*gg);
grossExtensivePlus = xbarPlus*LambdaPrimebarPlus;
selectionPlus = sum(deviationVec(negdeviationRange).*(LambdaPrime(negdeviationRange)-LambdaPrimebarPlus).*gg);

CountryCell.("intensive_h" + num2str(ii)) = F0*intensivePlus+(1-F0)*intensiveMinus;
CountryCell.("grossExtensive_h" + num2str(ii)) = F0*grossExtensivePlus + (1-F0)*grossExtensiveMinus;
CountryCell.("selection_h" + num2str(ii)) = F0*selectionPlus+(1-F0)*selectionMinus;
CountryCell.("sumAll_h" + num2str(ii)) = CountryCell.("intensive_h" + num2str(ii))+CountryCell.("grossExtensive_h" + num2str(ii))+CountryCell.("selection_h" + num2str(ii));
CountryCell.("LambdaPrimebar_h" + num2str(ii)) = F0*abs(LambdaPrimebarPlus)+(1-F0)*abs(LambdaPrimebarMinus);
CountryCell.("absLambdaPrimebarPlus_h" + num2str(ii)) = abs(LambdaPrimebarPlus);
CountryCell.("absLambdaPrimebarMinus_h" + num2str(ii)) = abs(LambdaPrimebarMinus);


if switch_figure
    figure(1);
    subplot(2,2,1); plot(deviationVec,hazardVec); title('Hazard function');
    subplot(2,2,2); plot(deviationVec,densityVec); title('Density');
    subplot(2,2,3); plot(deviationVec,LambdaPrime); title('\Lambda ''');
end